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! July  1,  1974  - June  30,  1975 

INTRODUCTION 

The  Polytechnic's  DARPA-STO- Aerophysics  Institute  program  during 
the  subject  period  was  centered  on  the  improvement  of  long-range  image 
recognition,  i.e.,  optical  discrimination,  and  on  the  study  of  certain  laser 
pnenomena. 

In  the  optics  area,  it  was  desired  to  improve  the  image  recovery  of 
an  object  of  finite  size  from  its  diffraction-limited  image  in  order  to  en- 
hance the  remote  observation  of  high-altitude  objects,  such  as  satellites, 
perceived  through  the  fluctuating  atmosphere.  A new  algorithm  is  proposed 
for  computing  the  transfer  of  a band-limited  function. 

In  the  laser  area,  the  major  task  was  the  screening  of  atomic  vapors, 
particularly  metal  vapors,  for  new,  efficient  lasers  in  the  visible  and  ultra- 
violet spectral  ranges,  using  computer-modeling  and  experimental  methods. 
Several  promising  candidates  are  identified.  In  another  vein,  the  theory  of 
unstable  resonator  cavitiss  was  studied  in  relation  to  lasers,  with  a view 
toward  improving  the  analytic  methodology  as  a complement  to  numerical 
techniques,  and  toward  improving  design  by  improving  the  understanding  of 
the  resonance  phenomena  in  cavity  configurations. 

The  developments  in  these  areas  arc  summarized  herein. 

A minor  effort  in  this  program  period  was  devoted  to  completion  of 
the  final  reportage  on  the  experimental  simulation  of  stratified  ocean  wakes 
in  our  thermally- stratified  air  wind-tunnel  (see  Appendix  III).  The  actual 
> work  had  been  completed  in  the  prior  program  period. 

The  principal  participants  in  the  program  were  Profs.  A.  Papoulis 
and  L.  Felsen  (currently  Dean  of  Engineering),  and  Dr.  William  Walter. 
Prof.  M.H.  Bloom  served  as  coordinator  of  the  program. 
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IMAGE-RECOVERY  OF  AN  OBJECT  Ot  FINITE  SIZE 
FROM  ITS  DIFFRACTION-LIMITED  IMAGE 

A.  Papoulic. 
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Under  the  usual  assumptions,  the  image  q(x,  y)  of  an  object  p(x,  v) 
can  be  expressed  as  a convolution 
00  00 

q(x,  y)  = j j p(x-5,  y-r?)  h(?,  p)  d?dp 
-00  -00 

In  the  coherent  case,  p(x,  y)  and  q(x,y)  are  amplitudes  and  h(x,y)  is  the 
point  spread  of  the  system.  In  the  incoherent  case,  p(x,y)  and  q(x,y)  are 
intensities  (illuminance  functions)  and  h(  x,  y)  is  the  square  of  the  point 

spread. 

Due  to  the  diffraction  effect,  the  transform  modulation  transfer 
function 

H(u,  v)  = J J h(x,  y)  e"j(uX+Vy*  dx  dy 
-00  -00 

is  space  limited: 

2 2 2 

H(u,  v)  = 0 for  u + v >a 

resulting  in  a distortion  of  the  object  p(x,  y). 

An  important  problem  in  image  processing  is  the  reduction  of  this 
distortion  (deconvolution).  Known  methods  involve  the  introduction  of  an 
optical  filter  on  the  aperture  plane  (apodization).  In  these  methods  the 
high  frequency  components  of  the  object  are  lost. 

In  this  work,  a method  is  developed  for  the  recovery  of  p(x,y)  under 
the  assumption  that  the  object  is  of  finite  size.  The  recovery,  in  p.inciple, 
of  p(x,  y)  is  a consequence  of  the  analytic  properties  of  its  transform 
P(u,  v).  However,  known  techniques  of  analytic  continuation  are  not  practical. 
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The  approach  proposed  here  is  based  on  a simple  efficient  algorithm 
involving  only  discrete  Fourier  series  and  Fast  Fourier  transforms. 

The  algorithm  involves  an  iteration  converging  to  the  exact  value  of  the 
object  p(x,  y)  in  the  absence  of  noise.  The  effects  of  noise  are  determined 
and  an  early  termination  of  the  iteration  is  suggested  for  optimum  results. 

Details  of  the  method  for  the  one-dimensional  problem  are  shown  in 


the  attached  paper. 


4 


IEEE  TRANSACTIONS  ON  CIRCUITS  AND  SYSTEMS,  VOL.  CAS-22,  NO.  9,  SEPTEMBER  i 975 

A New  Algorithm  in  Spectral  Analysis  and 
Band-Limited  Extrapolation 

ATHANASIOS  PAPOULIS,  fellow,  ieee 


Abstract — If  only  a segment  of  a function /(f)  is  given,  then  its  Fourier 
spectrum  F(to)  is  estimated  either  as  the  transform  of  the  product  of /(/ ) 
with  a time-iimitfd  window  wit) , or  by  certain  techniques  based  on 
various  a priori  assumptions.  In  the  following,  a new  algorithm  is  prop- 
osed for  computing  the  transform  of  a band-iimited  function.  The  algo- 
rithm is  a simple  iteration  involving  only  the  fast  Fourier  transform 
(FFT).  The  effect  of  noise  and  the  error  due  to  aliasing  are  determined 
and  it  is  shown  that  they  can  be  controlled  by  early  termination  of  the 
iteration. 

The  proposed  method  can  also  be  used  to  extrapolate  band-limited 
functions. 

I.  Introduction 

A CENTRAL  PROBLEM  in  Fourier  analysis  and 
spectral  estimation  is  the  determination  of  the 
tranefi  rm 

Fiw)  = f*  /(f)e'M  dt  (1) 

* - ® 

of  a signal  f(t)  in  terms  of  a finite  segment 

9(t)=mPT«),  Pi(t)  * |q’  |J|  > \ (2) 

of  f(t).  The  known  estimation  methods  can  be  separated 
into  the  following  two  classes. 

1)  (See  Blackman-Tukey  [1].)  The  available  segment 
g{t)  is  multiplied  by  a “window"  wit)  and  the  integral 

Fwi co)  = J w(r)/(Oe-'“'  dt 

= ir  Fid)  - y)Wiy)dy  (3) 

2rr  J_x 

is  used  as  the  estimator  of  Fiw).  The  pair  wit)  *-*  B'(cu) 
is  so  chosen  as  to  minimize  in  some  sense  the  resulting 
error  Fwiw)  - Fiw)  [2],  [3].  In  this  approach,  the  inverse 
transform/(/)H-(f ) of  the  estimator  Ffw)  is  zero  for  |f  | > T. 

2)  The  unknown  part  of  fit)  is  not  assumed  to  be  zero 
as  in  1)  but  it  is  determined  by  some  form  of  extrapolation. 
This  is  possible  only  if  various  a priori  assumptions  about 
for  form  of  Fiw)  are  made. 

Manuscript  received  November  15,  1974;  revised  March  26,  1 975. 
This  work  was  supported  in  part  by  the  Joint  Services  Electronics 
Program  under  Contract  F44620-69-C-0047  and  in  part  by  the  Ad- 
vanced Research  Projects  Agency  of  the  Department  of  Defense,  and 
was  monitored  by  ti^  Office  of  Naval  Research  under  Contract 
NOOOi  4-67-A-O438-O0i  7. 

The  author  is  with  the  Polytechnic  institute  of  New  York,  Farming- 
daie,  N.Y.  i i 735. 


Class  2)  includes  the  method  of  maximum  entropy  [4] 
(so.  Appendix  A)  and  the  te:hniques  for  extrapolating 
bat.d-limited  functions. 

Band- Limited  Functions 

We  shall  say  that  a function  fit)  is  band-limited,  if  its 
energy  E is  finite  and  its  transform  Fiw)  is  zero  outside  a 
finite  interval 

Fiw)  = 0,  for  M > a.  (4) 

A band-limited  function  is  analytic  in  the  entire  r-axis  [5], 
hence  its  Taylor-series  expansion  yields  fit)  for  every  t. 
We  hive  thus  an  extrapolation  method  (analytic  continua- 
tion) based  on  (4),  however,  as  it  is  well  known,  we  cannot 
use  it  for  rccovering/(f)  from  real  data. 

In  this  paper,  we  develop  a new  method  for  computing 
fit)  and  its  transform  Fiw)  from  the  given  segment  git). 
The  method  is  based  on  a simple  efficient  algorithm  involv- 
ing only  discrete  Fourier  series  and  a fast  Fourier  transform 
(FFT).  The  algorithm  is  presented  in  analog  form.  The 
numerical  implementation  is  carried  out  digitally  but,  since 
the  analog  to  digital  conversion  is  well  known,  the  details 
arc  omitted. 

We  note  that  the  proposed  algorithm  (Section  II)  has 
been  used  by  Gcrchbcrg  in  a recent  paper  on  picture 
processing  [6]. 

Another  method  of  extrapolating  band-limited  functions 
is  based  on  the  expansion  of  git)  into  a series  of  prolate 
spheroidal  functions  [7],  This  method  is  compared  with  ours 
in  Section  V. 

II.  An  Iteration  Method  for  Computing  the  Fourier 
Transform  of  a Bknd-Limited  Function 

We  wish  to  determine  the  Fourier  transform  Fiw)  of  a 
band-limited  function  fit)  (Fig.  1)  in  terms  of  the  finite 
segment  git)  of  fit)  [see  (2)].  The  proposed  method  is  an 
iteration  starting  with  the  transform 

Giw)  = G0iw)  = |r  git)e~Jut  dt 

of  the  given  segment  git)  = g0i 0-  The  nth  iteration  step 
proceeds  as  follows. 

We  form  the  function 

F„iw)  = G„.liw)p,iw),  p.iw)  «*  |q’  g > °g  (5) 


Copyright  © 1975  by  The  Institute  of  Electrical  and  Electronics  Engineer  Inc. 
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Fig.  1.  Unknown  pair/(r)  — F(w>  and  signals  of  first  iteration. 


obtained  by  truncating  G.-^cu)  as  in  Fig.  2,  and  compute 
its  inverse  transform 


/.(f)  = J do.  (6) 

We,  next,  form  the  function 

5.(0  =/»(*)  + [/(O  ~/.(0]Pr(0  ~ |?j  > T 

(7) 


obtained  by  replacing  the  segment  of /.(f)  in  the  interval 
(~T,T)  by  the  known  segment  g(t)  of /(f).  The  nth  step 
ends  by  computing  the  transform 

G.(o>)  = P g,(t)e~ta,dt  (8) 

* - OD 


of  the  function  g„(t)  so  formed. 

We  note  that  the  function  f„(t)  is  band-limited  and  it  is 
given  by 


m 


9,- 1(0 


sin  <xf 
rtf 


(9) 


This  follows  from  (5)  and  the  convolution  theorem  [8]. 
From  (9)  we  conclude  that/,(f ) can  be  obtained  as  the  output 
of  an  ideal  low-pass  filter  with  input  5„_,(f).  A two-input 
gate  at  the  output  of  the  filter,  switching  from  /.(f)  to 


Filltr 


MO 


Goo 


|»(0 


4.10 


Fig.  J.  Analog  simulation  of  iteration 


*4.10 


Fig.  4.  Computer  output  of  g,(i)  and  (7,(u>)  for /(/)  * sin  ai'ni»n<\ 
n - 2,  4,  6,  and  8. 


the  given  segment  y(f)  of /(f)  as  in  (7),  yields  5,(f)  (Fig.  3). 

It  will  be  shown  in  the  next  section  that  the  functions 
/(f)  and  F.( a)  tend  to  /(f)  and  F((o),  respectively,  as 
n - oo.  The  iteration  is  thus  a method  for  computing  F( d) 
and  for  extrapolating  a band-limited  function /(f). 

Numerical  Illustration 
We  have  applied  the  method  to  the  pair 

/(f)  = - — --  «-*  F(co)  - p,((o) 
rtf 

choosing  for  T the  value  rr/5cr.  The  numerical  results  are 
shown  in  Fig.  4 for  n = 2,4,6,8.  As  we  see  from  the  figure, 
although  the  data  interval  is  only  1/5  of  the  main  lobe  of 
/(f),  the  convergence  of  FJco)  to  the  unknown  spectrum 
p,(io)  is  remarkably  rapid.  With  this  T , the  method  ot 
windows  would  yield  no  satisfactory  estimate  of  F( cu).  We 
should  point  out,  however,  that  the  speed  of  convergence 
observed  for  this  example  is  due  to  the  fact  that  in  the 
expansion  (11)  of  the  function  sin  df/rrf,  only  the  first  few 
terms  are  significant.  As  we  show  in  Section  IV,  if  this  is 
not  the  case,  the  convergence  is  slow. 
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Mean-Square  Error  Reduction 
From  (5),  (7),  and  Parseval  formula  [$]  it  follows  that 


1 f* 
2ji  J 


|F(oj)  - G„_,(co)|J  do 


a: 

i: 
hi: 


I F(to)  - F.(co)|J  do 


1/(0  - m\ 


d,>S'_ 


I/O)  - ^(t)|J  dt 


I F(o)  - G,(oj)|j  do. 


(10) 


Thus  the  use  of  the  available  information  about  /(/) 
(band-limited  for  |co|  > a and  known  for  |r|  < 7”)  results 
in  a reduction  of  the  mean-square  error  twice  at  each 
iteration  step. 


III.  Convergence 

To  prove  the  convergence  of  the  iteration,  we  expand  the 
band-limited  function  f(t)  into  a scries 

fU)  = £ ak<f>k(t)  (11) 

k*-Q 

where  <pk(t)  are  the  prolate  spheroidal  wave  functions 
[see  (68)]. 

Theorem  1 

The  function  fm(t)  of  the  nth  iteration  is  given  by 

/.(')  =/(0  - f ak(l  - W*(0  02) 

k-0 


Applying  (13)  to  each  term  in  (11),  we  conclude  that 

/,(')  = I at[l  - (1  - WW 0 (17) 

k-0 

and  (12)  results. 

Theorem  2 
For  any  t, 

f,(t)  -*/(0  as  n - co.  (18) 

Proof:  As  we  see  from  (12),  the  error 

«.(0  = /(f)  -fft)  (19) 

is  given  by 

*„(0  = f a*(l  - W)  (20) 

k-0 

and  its  mean-square  value  En  by 

£.=  f"  f,J(0di=  i V(,-^)J"  (21) 

k-0 

[see  (70)].  We  maintain  that 

£„  -♦  0 with  n -*  co.  (22) 

Indeed,  the  energy 

E = i af  (23) 

k-0 

of /(r)  is  finite  by  assumption;  hence,  given  £ > 0,  we  can 
find  an  integer  N such  that 

E °k2  < e- 

t>.v 

The  eigenvalues  Xk  are  less  than  1 and  tend  to  zero 
monotonically  as  k -*  oo.  Therefore,  1 - Xk  < 1 and 


where  Xk  are  the  eigenvalues  of  (63)  corresponding  to  the 
eigenfunctions  <f>k(t). 

Proof:  Suppose,  first,  that 

AO  = <t>k(t). 

We  shall  show  by  induction  that 

/»(0  = 4A(0.  where  A„  = 1 - (1  - (13) 

Indeed,  (13)  is  true  for  n = I.  Suppose  that  it  is  true  for 
some  n ^ 1.  It,  then,  follows  from  (7),  (9),  and  (63)  that 

9 A)  = AMO  + (1  - A„)4>k(t)pT(t)-  (H) 

Inserting  into  (9)  and  using  (63)  and  (67),  we  obtain 

/.♦i(0  = ^k(t)  + (1  - A„)Xk<f>k(t)  (15) 

hence, 

fn+i(t)  = An+l  4>k(t) 

where 

^n+t  • Am  + (l  - AJ/.k,  n £ 1.  (16) 

Solving  the  recursion  equation  (16)  with  the  initial  condition 
A i = Xk,  we  obtain  (13). 


1 - Xk  Z 1 - ;,y,  for  k £ N. 

From  the  above  we  conclude  that 

£,  = E + E Vd  - h)1" 

k-0  *>.V 

< (1  - A.v)J"  £ V + Z V < (1  - Av)J"£  + E 

k-0  k> N 

and  (22)  results  because  £ is  arbitrary  and 

(1  - XN )2"  -»  0 with  n — oo. 

To  complete  the  proof  of  the  theorem,  we  observe  that 

\en(t)\  <,  (24) 

This  follows  from  the  fact  that  the  transform  of  e,(t)  is  zero 
for  |co|  > a [9],  Thus  e„(t)  -♦  0 with  n — x and  (18) 
results. 

We  note  for  later  use  that  [see  (7),  (12),  and  (66)] 

» hnT  T 

G(o)  = B E 4>k(bw),  B = — — , b = - 

k-0  \ <J  ( 7 

(25) 


(26) 
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FJL<0) 


Bp.( CO)  i 4s*  [1  - (1  - AJ%(M 


G,(cu)  = £.(o>)  + B £ atv^(l  - AJtykfM-  (27) 


fc-0 


Reasoning  similarly,  we  conclude  that  the  function  F„(co) 
and  G.(co)  tend  to  £(co)  with  n — co. 


IV.  Error  Analysis 

The  proposed  method  is  subject  to  the  following  errors: 
truncation  of  the  iteration;  aliasing  due  to  high  frequencies 
in  £(co);  noise;  computer  roundoff;  analog  to  digital 
conversion.  These  errors  depend  on  the  nature  of /(f)  and 
on  the  product  cT.  In  the  following,  we  examine  their  effect 
on  the  estimate  £.(cu)  of  £(e>).  The  analog  to  digital  con- 
version error  involving  the  evaluation  of  the  Fourier  in- 
tegrals by  discrete  Fourier  series  is  well  known  and  will  not 
be  discussed. 


Fig.  5. 


Firsi  15  eigenvalues  A.  of  the  integral  equaiion 
. sin  o(t  - r) 


*0- 


*(f  - t) 


dr  «=  XM.I). 


Truncation  Error 

If  the  iteration  is  terminated  at  the  nth  step  and  the  func- 
tion £.(co)  is  used  as  the  estimate  of  £(cu),  the  resulting 
mean-square  error  is  given  by  [see  (21)] 


*F(„)  kLIft.) 

1 

|M(«) 

/\ 1 A A/^ 

/ Y\  \ 

/ \ \ 
J y 

-A- 

“Tr — S f -*  * » » 

-cr  i 

D 9 

if' 

2*J-. 


Fig.  6.  Low-  anJ  high-frequency  components  of  F(w). 


| £(cu)  - £.(cu)|2  da)  = £ ak2(l  - Ak) 

*«a 


2ft 


(28) 


where  Ak  are  the  eigenvalues  of  (63).  As  we  see  from  Fig.  5, 
the  At’s  increase  as  the  product  oT  increases  and  they  tend 
to  1 with  oT  - x [10].  Thus  if  oT  is  large,  the  required 
number  of  iterations  is  small. 

The  error  £.  depends  also  on  the  coefficients  a..  If  n is 
small,  then  the  contribution  of  the  high-order  terms  in 
(28)  is  small.  However,  as  n increases,  their  effect  becomes 
relatively  more  significant  because  the  difference  1 - At 
tends  to  1 with  k -»  co. 

From  (27)  it  follows  after  some  thought  that  the  energy 
/.  of  G.(cu)  outside  the  band  (-a, a)  is  given  by 


must  be  considered  in  the  selection  of  the  cutoff  frequency  a. 
If  a is  large,  then  1 - A*  is  small;  however,  this  does  not 
lead  to  a reduction  of  the  error  £.  because  an  increase  of  a 
beyond  co0  results  in  an  increase  of  the  coefficients  ak.  If 
a is  smaller  than  cu0,  then  as  we  show  next,  aliasing  errors 
result.  The  optimum  it  is  a number  close  to  a)0  from  above. 


/,  = if  /G.(cd)|2  da  = t a»2Ak(  1 - A*)2"*1. 

2tt  J | o>i  >» 

Hence  [see  (21)], 


(29) 


k-0 


Aliasing  Error 

The  convergence  of /.(f)  to  the  function /(f)  is  based  on 
the  assumption  that  F(a)  = 0 for  |cu|  > a.  If  /(f)  is  not 
band' limited,  then  this  is  no  longer  true.  However,  as  we 
shall  show,  the  iteration  properly  terminated  yields  a useful 
estimate  of  £( co)  as  long  as  the  energy  £*  of  the  high- 
frequency  component  h(t)  of /(f)  is  sufficiently  small. 

We  introduce  the  functions  (Fig.  6) 


L{o>) 


/,  * \ i «k2(i  - A*)2"  = ^/4 

4 fc-o 


(30) 


}£(o)). 

|co|  ^ ff 

\o. 

|a>|  > a 

(°. 

|cu|  ^ a 

\Fico), 

|cu|  > cr 

(32) 


because  At(l  - At)  < i.  If,  at  the  nth-iteration  step,  the 
energy  /.  computed  from  (29),  is  large,  then,  the  energy  £. 
of  the  error  is  also  large.  However,  if  /,  is  small,  £,  is  not 
necessarily  small.  The  inequality 

- . A. 


+ I flk2 

k>n 


(31) 


1 - A. 

obtained  readily  from  (21)  and  (29),  gives  an  upper  bound 
for  £.+ , but  this  is  useful  only  if  some  information  about  the 
speed  of  convergence  of  the  series  (11)  is  available. 

Suppose,  finally,  that  £(co)  = 0 for  |a>|  > a>0  ^ut  l^e 
precise  value  of  o0  is  not  given.  In  this  case,  the  following 


H(  a)  = 

and  their  inverse  transforms  /(f)  and  h(t).  Clearly, 

£(  oi)  = I(o>)  + H(d) 

fit)  = lit ) + MO- 

The  iteration  scheme  is  obviously  linear  in /(f);  therefore, 
/.(f)  can  be  written  as  a sum 

m - /„(o  + Kit)  (33) 

where  /.(f)  and  /i.(f)  are  the  contributions  due  to  /(f) 
and  MO.  respectively.  As  we  see  from  (32),  the  function 
/(f)  is  band-limited;  hence,  the  component  /.(f)  of  /.(f) 
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tends  to  /(f)  with  n -*  oo.  It  suffices,  therefore,  to  examine 
the  behavior  of 

We  expand  the  function  h(t)  into  a series  in  the  interval 

(~T,T) 

h(t)  = i c„M f),  Ul  ^ T <34> 

k-0 

[see  (69)].  Clearly, 

i ck2xk  = f h2(t)dt  < £* 

k-0  J -T 

= -1  I*  |F(co)|2  dco.  (35) 
2it  J |w| >» 


It  is  not  difficult  to  see  from  (75)  that 


E{dkd*}  = 

|S/Ak,  k - r 

(0,  k ft  r 

(43) 

£{|*,(0))|2}  = 

2Sitb  £ ^*k2(6w)  = 2ST 

k-0 

(44) 

The  nth-iteration  function  X,(co)  is  given  by  the  nght 
side  of  (26)  where  the  coefficients  ak  are  replaced  by  <f, 
Hence,  its  mean-square  value  R , is  given  by 

Rn  = £{|X.(o,)|2} 

= 2Snb  £ f1  ~ ° ~ — 1 <>k3iboi) 
k-0  l J 


From  (17)  it  follows  that 

/!„( t)  = f ck(l  - pk')U'),  Mk  = 1 ~ <36) 

k-0 

The  function  hn(t)  is  band-limited  by  construction  and  its 
energy  equals 


£*  = i ^k2d  - Mk-)2 

k-0 

= £ Ck2A*2(l  4-  Hk  + •••  + ^r1)2-  (3?) 

k-0 


The 

numbers  Ak  and  pk  are  positive 

and  less  than  1 ; 

hence. 

1 + pk  + • • • 4-  i4~ 1 < 

n 

and 

£«k  < "2  I ck2A»  = "2 

k-0 

Ek.  (38) 

The  above  bound  is  too  pessimistic;  the  actual  value  of 
might  be  considerably  smaller. 


None  and  Roundoff  Error 

Suppose  that  for  |f|  < T,  the  known  signal  is  a sum 

g(t)  + y(  0 (39> 

where  y(t)  is  a random  process  with  autocorrelation 
SW'iW'j)}  = Wuh)- 

At  the  nth  iteration,  f„(t)  is  replaced  by/„(f)  + *„(0  where 
the  random  component  xR(f)  is  due  to  j>(r).  We  shall 
assume  that  the  noise  y(t)  is  white  as  in  (78) 

£(f„f2)  = Spr(0<5('i  ~ h)-  ^ 

We  expand  the  transform 

Y{(0)  = |T  y(t)e-Ja ' dt  (41) 

of  j>(f)  into  a Karhunen-Loeve  series  in  the  interval  (-0,0) 
[see  (74),  (80),  and  (25)] 

XJw)  - = Bpt((o)  £ dk\l Xk  <f>k(bco).  (42) 

k-0 


= 2Snb  £ (1  + Pk  + + ffk~')2<t> k2(M-  (45) 

k-0 

This  error  is  similar  to  the  aliasing  error  £„»  in  (37)  and  it 
is  bounded  by 

Rn  <;  2STn}  (4«) 

as  in  (38). 

The  preceding  analysis  can  be  used  to  determine  the 
effect  of  the  roundoff  error  on  the  estimate  F,(a>)  of  F((o). 

In  the  evaluation  of  the  Fourier  transform  G(w)  of  g(t) 
by  a sum,  roundoff  errors  are  introduced.  If  these  errors 
are  independent  with  standard  deviation  a (smallest  available 
digit),  then,  their  effect  on  £„(co)  is  approximately  equivalent 
to  the  presence  of  a noise  term  y(t)  as  in  (39)  with  .5  — 
a2S  where  5 is  the  sampling  interval.  The  resulting  mean- 
square  error  at  the  nth  iteration  is  given  by  R „ as  in  (45). 
Similarly,  the  roundoff  error  introduced  at  the  fcth-iteration 
step  causes  a mean-square  error  in  the  evaluation  of  £„(a>) 
equal  to  R„-k- 

V.  Comparison  with  the  Frol  ate  Spheroidal  Expansion 

Method 

We  shall  discuss  briefly  the  method  of  prolate  expansion 
for  dr’ermining  E(cu)  from  g(t)  and  will  compare  it  with  our 
approach. 

We  expand  the  tir  ♦-limited  function  g(t)  into  a senes 
as  in  (69).  The  coefficients  ak  are  given  by 

a*  = f siOUO  dt.  (47) 

Ak  J -T 

We  insert  the  numbers  so  obtained  in  the  expansion  (11). 
The  resulting  sum /(f)  equals  the  bard-limited  extrapolation 
/(f)  of  the  piece  g(t).  The  transform  F(co)  of /(f)  is  obtained 

from  (25) 

£(o>)  = B £ M £ * (*8) 

° VAk 

Thus  to  evaluate  £(co)  with  this  method,  we  musv  store 
the  sample  values  of  the  functions  <fik{t)  and  use  them  to 
evaluate  the  integral  in  (47).  Our  method  involves  only 
discrete  Fourier  series  and  makes  use  of  available  algorithms 
and  the  N log  N economy  of  the  FFT. 
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In  most  applications,  the  product  oT  is  large.  This 
assumption  does  not  reduce  the  computations  in  the  prolate 
method.  In  our  method  it  leads  to  a rapid  convergence 
that  can  be  sensed  by  computing  E 
If  noise  or  aliasing  errors  are  present,  no  method  permits 
the  precise  evaluation  of  F(oi).  An  approximate  result  is 
possible  with  the  prolate  method  if  the  series  (48)  is  trun- 
cated. In  our  approach,  the  iteration  must  be  carried  out  to  a 
value  n such  that  the  errors  Enh  and  Rn  in  (37)  or  (45)  are 
sufficiently  small.  If  n is  too  small,  the  truncation  error  £„ 
[see  (28)]  might  be  unacceptable.  In  a specific  example,  a 
compromise  must  be  sought;  the  bounds  (JO),  (38),  and  (46) 
give  some  guidance. 

VI.  Relationship  to  Frequency-Dependent  Windows 
and  Integral  Equations 

In  the  proposed  iteration,  the  function  Ffai)  of  the  /ith 
step  is  a linear  functional  of  the  segment  g(t ^ We  denote  by 
y(a), t)  the  limit  of  F„{(o)  corresponding  to  g(t)  = <5(f  - t). 
Since 


Taking  Fourier  transforms  of  both  sides  of  (53)  we  obtain 
(49)  where 

y(co,x)  = I*  k(t,x)e~iM  dt.  (54) 

J - ao 

The  discrete  form  of  (52)  is  a system  of  linear  equations. 
Appendix  A 

The  Method  of  Maximum  Entropy 
In  the  numerical  evaluation  of  F(co),  the  integral  in  (1)  is 
replaced  by  a sum 

F((o)  = £ /.  = «/("«)■  (55) 

■ ■-  a. 

This  is  the  exact  if  F( to)  = 0 for  |tu|  > n/a.  Otherwise,  it  is 
approximate  subject  to  aliasing  errors.  Since  only  the 
segment  g(t)  of  /(f)  is  known,  the  above  sum  must  be 
estimated  in  terms  of  its  2 M + 1 coefficients 

/„  = = ag(na),  |n|  £ M = [T/a].  (56) 


9(0  = J 9(0&(t  “ 0 dx 

it  follows  from  the  linearity  of  the  iteration  that 
F(w)  = | g(x)y(co,x)  dx. 


With 


(49)  yields 


w(ci),t)  = yi(t),t)eJa 


F((o)  = JT  M.oi,t)g(t)e-ia“  dt. 


(49) 


(50) 


This  is  of  the  form  (3)  with  the  difference  that  the  window 
w(a),f)  is  frequency  dependent.  Unlike  the  estimate  Fw((d) 
in  (3),  the  inverse  of  F(wi)  in  (50)  is  not  zero  for  |f|  > T 
but  it  equals  the  unknown  segment  of /(f). 

We  shall  next  show  that  the  function  >( cu.r)  in  (49)  is  the 
Fourier  transform  of  the  resolving  kernel  of  an  integral 
equation.  Since  /(f)  = g(t)  for  |f|  < T,  the  Fourier 
trar.  iform  F(co)  of /(f)  can  be  written  as  a sum 


Ffjto)  — j*  g(t)e-ia' dt  + \ /(Oe-^dt.  (51) 

J|f|<r  *M>l>r 


The  method  of  maximum  entropy  is  developed  for  the 
estimation  of  power  spectra  [11],  [12].  In  this  case,  F(to) 
is  real  and  positive.  The  function  F(to)  is  so  selected  as  to 
maximize  the  integral  (entropy) 

log  F(o>)  dto 

J -*!• 

subject  to  the  2M  + 1 constraints  (56).  It  is  easy  to  show 
that  [13],  this  assumption  leads  to  the  conclusion  that 
F(to)  must  be  rational  with  no  finite  zeros.  Since  F(oi)  £ 0, 
it  follows  from  the  Fejer-Ricsz  theorem  [5]  that 


F(co)  = 


M 


(57) 


1 + I 


-jUa\ 


k - 1 


where  all  the  zeros  z,  of  the  polynomial 


M 


Y(0  = i + I n*' 

k-  1 


are  inside  the  unit  circle. 

Our  problem  is  thus  reduced  to  the  determination  of  the 
M + 1 constants  p,ylt-  ■ -,ym.  The  number  of  independ;nt 
equations  in  (56)  is  also  M + 1 because/.,  =/,*• 

We  maintain  that 


But  F((o)  = 0 for  |a>|  > a by  assumption;  hence, 


for  |cu|  > <r.  (52) 

This  is  a linear  integral  equation  for  the  unknown  seg- 
ment of  /(f).  Denoting  by  k(t,x)  its  resolving  kernel,  we 
obtain 

/(<)  = I*  9(0H f,x)  dx.  (53) 


I S.-kFk  = -9»  n = l,  - M (58) 

k-  1 

and 

P - g0  + I 0-kFk-  (59> 

k-1 

Thus  to  find  the  yk's  we  solve  the  system  (58).  The  special 
form  gH. k of  the  coefficient  matrix  (Toeplitz)  simplifies  the 
computations.  Inserting  the  solution  into  (59),  we  obtain 
the  constant  p. 


PAP0UU3'.  ALGORITHM  IN  SPECTRAL  ANALYSIS 


10 


Proof:  With 


"(*)-  I U 


wc  have  [see  (55)]  F(co)  = H(e Ja‘)  and  (57)  yields 


H(z)  (l  + 1 


P 

i + I y*zk 

*-  I 


The  zeros  1/z,*  of  the  denominator  F*(l/z*)  of  (61)  are 
outside  the  unit  circle,  hence,  the  inverse  z-transform  of  the 
right  side  of  (61)  is  zero  for  n > 0 and  it  equals  p for  n = 0. 
Equating  the  inverse  z-transforms  of  both  sides,  we  obtain 

/,  + = (o,’  n > O'  (62) 

The  above  yields  (58)  for  n = 1 • - ,M  [see  (56)]  and  (59) 
for  n = 0. 

We  note  that  the  maximum  entropy  method  can  be 
considered  as  a method  of  estimating  F(oj)  in  terms  of  gn 
with  the  a priori  assumption  (Class  2))  that  F( co)  is  of  the 
form  (57). 

Appendix  B 

The  Prolate  Spheroidal  Functions 
The  eigenfunctions  <f>k(t)  of  the  equati  jn 

fr  m 

J -t  n(t  - t) 

known  as  prolate  spheroidal  wave  functions  [7],  are 
orthogonal  in  the  intervals  (-00,00)  and  (~T,T) 


HOMO  dt  = 


1,  i 


\T t uouo  dt  = {£• 


i = fc 
i * fc' 


A band-limited  function  /(f)  can  be  expanded  into  a 
series  in  the  interval  (-00,00) 

f(0  = f a WO 

k*0 


r co 

Ok  - I /(two  dt. 

J - 00 


An  arbitrary  function  g{t ) can  be  expanded  into  a series  in 
the  interval  (-  7*,7*) 

9(0  = f ak4>k(t) 

*- 0 


7 r 9(tW0  dt. 

*k  J -T 


From  the  above  and  (64)  we  conclude  that 


f*  00  00 

f2[t)dt=  l a,1 
J -00  *-0 

f g2(t)  dt  - f a*’A, 

J - J-  km0 


The  eigenvalues  /*  are  such  that 
1 > A0  > • • • > A*  > • • • > 0 and  A*  ->  0, 

with  k -*  00.  (65) 

The  functions  <f>k(t)  equal  their  own  Fourier  transforms 
properly  truncated  and  scaled: 

B 

4>k(0  ~ 77  <Pk(b(o)p,(oj), 

B = \ 2nTlo,  b = T/o 
4>k(t)PT(t)  ^vAt  (f>k(bw).  (66) 

Thus  the  solutions  of  (63)  are  band-iimited  hence, 


We,  finally,  note  that 

f 4>*\t)  dt  - a/n.  (71) 

k-0 

This  follows  by  expanding  the  band-limited  function 
sin  o(t  - x)/n(t  - x)  into  a series  as  in  (68)  and  using  (63). 

Appendix  C 

The  Karhunen-Lc  eve  Expansion 
Consider  a stochastic  process  y(t)  with  autocorrelation 

RUrdi)  = £{.K', )>’*('*)}• 

With 

r(u,t>)  = £{F(u)r*(u)} 

the  autocorrelation  of  the  Fourier  transform  Y(co)  of  y{t), 
it  is  easy  to  show  that  [!4] 


P 00  poo 
" — oo  " — a 


R(tltt2)  exp  (-/«»,  - vt2))  dtx  dt2. 


With  r(u,t>)  as  above,  we  form  the  equation 


/>• 


v)\j/(v)  dv  = y \j/(u) 


and  denote  by  fj/k(u)  and  yk  its  eigenfunctions  and  eigen- 
values. It  can  be  shown  that  in  the  expansion 

T(a>)  = £ z*^*(cu),  H £ a (74) 

*-0 

the  random  variables  zk  are  orthogonal  [14]  and 


£{z,z**} 


f y*.  «'  « k 

10,  i # k' 
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Nonstationary  white  noise:  We  say  that  the  process  .p(f) 
is  white  noise  if 

RUi,t2)  = qUl)6(tl  - h)  (76) 

where  q(t)  £ 0 is  its  “average  intensity.”  With  Q(a>)  the 
Fourier  transform  of  q(t ),  it  follows  from  (72)  that 

r(u,u)  = f\u  - v).  (77) 

If  y(t)  = n(t)pT(t)  and  n(f)  is  stationary  white  noise  with 
autocorrelation  S5(f,  — t2),  then 

ROi.ti)  = SpT(tl)6(tl  - t2).  (78) 

In  this  case, 

q(t)  = Spr(t)  Q((o)  = 2 S sin  7ct>/£t>.  (79) 

From  the  above  and  (77)  it  follows  that  the  integral 
equation  (73)  takes  the  form 

2S  f sin  T(u  -_v)  ^ (80) 

u - v 

As  we  see  from  (63),  its  eigenfunctions  i/t^u)  equal  the  prolate 
spheroidal  functions  <t>k(t)  properly  scaled. 
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UNSTABLE  RESONATORS 


L.  Felsen 

As  noted  in  the  proposal,  we  have  developed  a new  method  of 
analysis  of  empty  unstable  resonator  cavities  whereby  the  finite  mirror 
configuration  is  regarded  as  a waveguide  in  the  direction  transverse  to 
the  resonator  axis.  Invoking  coupling,  at  the  mirror  edges,  between  two 
appropriate  waveguide  modes,  we  have  been  able  to  explain  in  a quantitative 
manner1  the  intricate  behavior  of  the  eigenmode  loss  curves,  which  have 

been  furnished  in  the  literature  by  numerical  solution  of  the  resonator 
2,  3 

integral  equation 


The  effort  during  the  current  period  has  been  directed  along  two 
lines:  a)  to  further  solidify  the  basic  approach  so  as  to  make  it  applicable 
also  to  mirror  configurations  other  than  the  hyperbolic  shape  considered 
so  far;  and  b)  to  extend  the  analysis  to  resonators  filled  with  inhomogeneous 

and  active  material. 

On  item  a),  we  have  completed  a detailed  study  of  the  modal 
resonance  equation  to  determine  the  eigenmode  loss  characteristics  of 
symmetrical  resonators  with  cylindrical  and  spherical  mirrors  of  hyper- 
bolic shape.  We  have  then  shown  that  the  waveguide  model  can  be  applied 
as  well  to  other  mirror  configurations  (asymmetrical,  telescopic,  etc.) 
and  that  the  properties  of  these  structures  can  be  related  to  those  of  the  . 
symmetrical  hyperbolic  type.  The  equivalence  relations  so  derived  from 
the  waveguide  model  agree  with  those  presented  in  the  literature  from 
considerations  involving  the  resonator  integral  equation.  This  serves  as 
further  substantiation  of  the  validity  of  the  waveguide  model.  The  results 
of  this  study  have  been  prepared  in  the  form  of  a manuscript  which  is  to  be 
published  in  the  Applied  Optics  Issue  of  the  Proceedings  of  the  1975  CLEA  . 

An  essential  part  of  the  waveguide  approach  is  the  separation  of  the 
overall  resonator  behavior  into  "canonical"  constituents  which  highlight  the 
separate  roles  played  by  edge  diffraction  and  by  propagation  between  the 
edges  and  the  resonator  interior.  The  geometric  optical  nature  of  the  edge 
diffraction  process,  and  its  influence  on  the  eigenmode  fields  and  losses, 
has  been  recognized  and  conjectured  but  has  not  previously  been  incor- 
porated into  a quantitative  theory.  We  have  applied  and  generalized  a pre- 
viously developed  ray-optical  theory  for  scattering  by  edge  discontinu 


in  waveguides  to  show  the  adequacy  of  the  geometric-optic  edge  diffraction 
mechanism  at  large  Fresnel  numbers.  We  have  also  shown  the  transition 
from  this  regime  to  that  of  small  Fresnel  numbers  where  edge  diffraction 
is  more  appropriately  described  by  a local  open-ended  parallel  plane  wave- 
guide model.  These  new  results  accommodate  the  entire  range  of  para- 
meters from  large  to  small  Fresnel  numbers  within  the  same  quantitative 
format  based  on  the  waveguide  approach. 

On  item  b),  work  has  been  initiated  on  two-dimensional  parallel 
plane  resonators  filled  with  a medium  that  is  inhomogeneous  in  the  direction 
transverse  to  the  resonator  axis.  The  parallel  plane  configuration  was 
chosen  initially  to  permit  separability  of  the  wave  equation  in  the  inhomo- 
geneous ly  filled  region.  A principal  goal  of  this  study  is  to  assess  the 
separate  influence  of  the  mirror  edges  and  the  medium  profile  (both  refrac- 
tive index  and  gain)  on  the  resonator  behavior.  The  problem  is  being  for- 
mulated for  arbitrary  (though  slowly  varying,  over  a wavelength  interval) 
profile  shapes  by  recourse  to  WKB  and  uniform  asymptotic  methods  avail- 
able from  the  theory  of  differential  equations. 


References 

1.  L.W.  Chen  and  L.  B.  Felsen,  IEEE  J.  Quant.  Electron.  QE-9,  1102,  1973. 

2.  R.L.  Sanderson  and  W.  Streifer,  Appl.  Opt.  2129  (1969). 

3.  A.E.  Siegmanand  H.Y.  Miller,  Appl.  Opt.  9,  2729  (1970). 

4.  C.  Santana  and  L.  B.  Felsen,  " Unstable  Resonators  - Two-Dimensional 

and  Three-Dimensional  Losses  by  a Waveguide  Analysis  , .o  be 
published  in  CLEA  Special  Issue  of  Applied  Optics. 

5 D.V.  Batorsky  and  L.  B.  Felsen,  Radio  Sci.  6.,  911  (1971). 


14 


r 


METAL  SCREENING  PROGRAM 

Vapors  of  atoms  with  suitable  low-lying  energy  level  structures  can 
provide  active  media  for  efficient,  pulsed  gas-discharge  or  possibly  chemical 
lasers.  1 The  goal  of  the  metal  screening  program  is  to  investigate  these 
metallic  vapors  as  the  active  media  for  more  efficient  lasers  in  the  visible 
and  ultraviolet  spectral  regions.  A three  phase  program  is  being  pursued  con- 
sisting of: 

1.  System  Evaluation 

2.  Computer  Modeling 

3.  Experimental  Investigation 

1.  In  the  initial  evaluation  the  criteria  developed  earlier  1 concerning 
the  ordering,  spacing  and  parity  of  energy  levels,  atomic  transition  probabili- 
ties and  branching  ratios  were  applied  to  all  the  elements.  The  most  promising 
possibilities  for  the  ultravic’et  and  blue-green  spectral  regions  are  listed  below: 


i 


Ultraviolet  Candidates 

2414  X Selenium 
2770  X Tellurium 

Blue-Green  Candidates 

4685  X Germanium 
4722  X Bismuth 
4942  X Chromium 

These  candidate  laser  transitions  were  selected  for  two  reasons:  first, 
their  wavelengths  are  in  spectral  regions  of  strong  interest;  and  second,  the 
2414  X selenium,  2770  X tellurium  and  4942  X chromium  lines  have  transition 
probabilities  which  are  smaller  than  that  of  the  similar  5105  X copper  vapor  laser 
transition.  A low  value  of  the  proposed  laser  line's  atomic  transition  proba- 
bility is  important  because  of  the  superradiance  limitation  on  this  class  of  pulsed 
self-terminating  or  cyclic  lasers.  For  the  cyclic  laser  in  copper  vapor  we  have 
measured  a small  signal  gain  coefficient  of  4 dB/cm  on  the  5105  X transition. 

With  such  a large  gain,  superradiance  (or  better  superfluorescence)  is  the  limit- 
ing process  in  scaling  the  copper  vapor  laser.  If  cyclic  lasers  can  be  developed 
on  transitions  with  smaller  atomic  transition  probabilities,  higher  densities  or 
larger  volumes  could  be  utilized  before  reaching  this  limit.  The  proposed  cyclic 
laser  transition  at  4942  X in  chromium,  for  example,  has  a transition  proba- 
bility of  4 x 105  which  is  five  times  lower  than  that  of  copper. 

^ 
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The  atomic  energy  levels  of  Moore  and  the  atomic  transition  proba- 

3 

bilities  of  Corliss  and  Bozman  have  been  used  in  this  initial  evaluation.  Al- 
though the  work  of  Corliss  and  Bozman  is  the  broadest  treatment  of  atomic 
transition  probabilities  presently  available,  errors  as  large  as  a factor  of  20 
have  been  found  when  their  values  have  been  compared  with  more-accurate 
measured  values.  Because  the  Corliss-Bozman  values  remain  the  only  con- 
sistent set  of  transition  probabilities  spanning  the  stronger  spectral  lines  of 
m Jst  of  the  solid  elements,  there  has  been  considerable  interest  in  attempting 
to  improve  these  values. 

Accurate  values  of  atomic  transition  probabilities  are  most  important 
in  the  evaluation  of  new,  efficient  laser  systems.  Therefore  as  part  of  the 
metal  screening  program  we  have  explored  the  possibility  of  improving  the 
Corliss-Bozman  values  of  atomic  transition  probabilities.  A procedure  was 
developed  to  obtain  improved  transition  probabilities  from  the  spectral  line 
intensities  of  Meggers,  Corliss  and  Scribner  by  applying  a correction  to  the 

temperature  of  the  Meggers-Corliss-Scribner  arc  and  by  incorporating  each 
element's  partition  function  in  the  procedure.  This  procedure  is  described 
in  Appendix  I.  In  Appendix  I the  temperature  of  the  Meggers-Corliss-Scribner 
arc  is  shown  to  be  closer  to  6200°K  than  the  5100°K  they  assumed.  Recently 

4 

the  NBS  group  published  a revised  set  of  spectral  line  intensities.  These  re- 
vised intensities  are  now  being  analysed  to  refine  the  procedure  and  develop 
a table  of  improved  atomic  transition  probabilities  for  transition  of  interest 
as  potential  efficient  laser  candidates. 

2.  A computer  model  of  the  dynamics  of  these  cyclic  lasers  has  been 
developed  and  is  now  being  refined.  The  goal  is  to  increase  our  understanding 
of  the  importance  of  various  physical  processes  involved  as  well  as  to  provide 
insight  on  scaling  these  lasers.  The  computer  model  is  described  in  Appendix 
II.  Thus  far  the  model  has  been  applied  to  two  metal  vapor  laser  systems: 

1)  the  copper  vapor  system  - to  compare  the  model's  predictions  with  the 
actual  performance  of  our  copper  vapor  lasers,  and  2)  the  bismuth  vapor  system- 
to  analyse  the  effect  of  andassessthe  importance  of  processes  which  compete 
with  the  basic  laser  dynamics. 

During  the  development  phase  of  the  computer  model,  the  model  was 
tested  on  two  copper  vapor  laser  geometries.  Qualitative  agreement  was  ob- 
tained between  the  model  predictions  and  the  experimental  measurements  but 


16 


a detailed  comparison  has  not  yet  been  carried  out.  Attention  was  shifted  to 
the  bismuth  vapor  system  where  cyclic  laser  action  has  not  been  observed  at 
4722  A in  spite  of  demonstration  attempts  at  various  laboratories  including 
the  Polytechnic.  Because  bismuth  has  three  valence  electrons,  there  are 
more  excitation  paths  available  in  a bismuth  discharge  than  in  the  single 
valence  electron  system  in  copper.  Cascade  filling  of  the  lower  laser  level 
is  also  more  likely  to  be  a detrimental  process  in  bismuth  vapor. 


To  examine  the  question  of  whether  these  competing  processes  could 

account  for  the  absence  of  laser  action  in  bismuth  vapor,  the  computer  model 

was  modified  to  include  electron  excitation  to  and  spontaneous  radiation  from 

2 

fifteen  levels  of  the  bismuth  atom.  Atomic  energy  levels  of  Moore  and 
transition  probabilities  of  Corliss  and  Bozman  were  used.  Electron  excitation 
and  ionization  cross  sections  were  estimated  using  the  semiclassical  formulas 
of  Gryzinski.  ^ 


For  our  experimental  conditions  the  model  predicted  that  a population 

inversion  would  be  achieved  and  laoer  action  at  4722A  would  result  in  the 

bismuth  vapor  discharge.  The  competing  processes  therefore  did  not  appear 

to  be  sufficiently  detrimental  to  prevent  laser  action.  In  May  of  1975  the 

7 

electron  impact  measurement  by  Williams,  Trajmar  and  Bozinis  (WTB) 
of  the  electron  excitation  cross  section  at  40  eV  of  the  upper  an{l 

lower  bismuth  levels  became  available.  At  40  eV  our  Gryzinski  esti- 

mate of  the  cross  section  of  the  upper  laser  level  (^Pj^)  was  aPP r oxi rn  at  e ly 

one-half  that  of  the  WTB  value  and  our  Gryzinski  estimate  of  the  cross  section 

2 

for  electron  excitation  to  the  lower  laser  level  ( D-/,)  was  approximately 

1 5/  t 

/ 30  that  of  the  WTB  value.  We  then  used  these  WTB  values  to  normalize 
our  relative  Gryzinski  excitation  curves;  that  is,  we  increased  our  upper  laser 
level  excitation  rates  by  a factor  of  2 and  our  lower  level  excitation  rates  by 
a factor  of  30  throughout  the  energy  range.  Still  the  computer  model  predicted 
laser  action  at  4722  A with  a 30  kW  peak  output  power  under  conditions  which 
correspond  to  thosecf  our  apparatus.  Thus  far,  then,  competing  processes  do 
not  appear  to  account  for  the  absence  of  laser  action  in  bismuth  vapor. 


f 
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The  computer  model  is  now  being  refined.  Reaction  processes,  in- 
cluding dimer  formation  and  multiple  ionization,  are  being  incorporated  into 
the  model.  Attention  is  also  being  directed  to  the  initial  nanosecond  during 
which  the  electron  temperature  in  the  model  rises  to  unphysically  high  values 
(106  eV).  The  interplay  of  experimental  results  with  computer  model  predictions 
greatly  facilitates  the  progress  in  understanding  the  importance  of  the  physical 
processes  involved. 

3.  Experimental  work  is  underway  on  both  the  selenium  and  the 
germanium  vapor  systems.  Selenium  is  strongly  electronegotive,  lying  under 
oxygen  in  the  periodic  table.  When  the  vapor  pressure  of  selenium  reached 
-lOmtorr  in  our  conventional,  longitudinal-discharge,  metal  vapor-rare  gas 
mixture  systems,  the  discharge  extinguished  and  could  not  be  reestablished 
with  voltages  as  high  as  20  kV.  Shorter  longitudinal-discharges  were  tried 
without  success.  Obviously  systems  with  higher  values  of  E/p  are  required. 

We  are  now  constructing  a transverse-discharge  system  to  overcome  this 
problem.  Shaped  electrodes  have  been  constructed  and  preionization  will  be 
utilized  to  establish  a uniform,  large-area  discharge.  The  initial  testing  of 
this  new  transverse-discharge  system  will  be  carried  out  at  room  temperature 
with  a nitrogen  fill.  The  pulsed  nitrogen  laser  at  3371  X is  a similar  self- 
terminating transition  which  can  be  used  to  test  and  calibrate  the  new  apparatus. 

Our  initial  examination  of  germanium  vapor  for  laser  action  at  4685  X 
has  been  quite  promising.  Unlike  bismuth,  germanium  is  transparent  at  the 
proposed  laser  wavelength.  Blocking  the  back  mirror  produces  a decrease  in 
the  detection  signal  at  4685  X.  The  hot  zone  available  in  the  metal  vapor  system 
used  was  - 12  cm  long.  The  germanium  vapor  system  will  next  be  examined 
in  a metal  vapor  system  where  a longer  hot  zone  is  available  (~  80  cm). 

Chromium  is  also  a promising  candidate  to  test  in  the  longer  metal 
vapor  facility.  As  alre.  , mentioned  above,  the  Corliss-Bozman  value  for 
the  transition  probability  of  the  4942  X proposed  laser  line  is  5 times  smaller 
than  that  of  the  green  5105X  transition  in  copper  vapor.  Therefore  the  supcr- 
radiance  limitation  will  be  less  severe  in  chromium  than  in  copper. 
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In  conclusion  then,  a three  phase  effort  is  being  carried  out  in  the 
screening  of  metal  vapors  for  efficient  visible  and  ultraviolet  laser  systems. 
System  evaluation  has  identified  several  promising  candidates.  A procedure 
has  been  developed  to  improve  the  values  of  atomic  transition  probabilities 
determined  from  spectral  line  intensities.  A computer  model  of  these  pulsed 
metal  vapor  systems  has  beon  developed.  The  model  is  being  refined  to  de- 
termine the  importance  of  various  physical  processes,  complete  the  evaluation 
of  the  candidate  metallic  vapor  systems  and  examine  the  scaling  of  these 
systems.  Experimentally  selenium  and  germanium  vapor  systems  are  under- 
going further  examination  - selenium  in  a transverse  discharge  and  germanium 
in  a longer  longitudinal  discharge.  Chromium  will  also  be  examined  in  the 
longer  longitudinal  discharge.  It  is  recommended  that  this  work  be  supported 
for  another  year  to  carry  out  these  evaluations. 
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IMPROVED  VALUES  OF  ATOMIC  TRANSITION  PROBABILITIES 
W.  T.  Walter  and  N.  Solimenc 

The  rate  at  which  an  appropriately  excited  atom  or  molecule  spontaneously  radi 
ates  in  a particular  spectral  line  is  characterized  by  a transition  probability, 

Aujf(Sec  w,lich  is  sometimes  referred  to  as  the  Einstein  A coefficient.  The  wave- 
length of  the  spectral  line  is  related  to  the  energies  of  the  upper  (u)  and  lower  (f ) levels, 
X = he/  (E^  - Ej).  Accurate  values  of  these  transition  probabilities  arc  required  when- 
ever the  kinetics  of  processes  involving  excited  states  of  atoms  or  molecules  are  ex- 
amined; for  example,  in  the  search  for  more  efficient  visible  and  ultraviolet  lasers,  in 
the  study  of  chemical  reactions,  for  astrophysical  determinations  of  solar  abundances, 
interstellar  cloud  composition,  etc. 

Transition  probabilities  can  be  both  calculated  and  measured.  They  can  be  calcu- 
lated from  first  principles  only  for  the  lightest  elements,  but  the  calculations  may  be 
extended  to  heavier  elements  by  including  experimental  parameters.  Experimental 
measurements  of  transition  probabilities  can  be  divided  into  two  classes;  (1)  methods 
which  require  knowledge  of  the  density  of  the  radiating  species,  such  as  the  hook 
method  of  Rozhdestvenskii  and  emission  and  absorption  techniques,  and  (2)  methods 
which  are  independent  of  the  vapor  density  such  as  level  crossing,  optical  double  reso- 
nance and  phase  shift  techniques.  Only  these  latter  measurement  methods  can  be  relied 
upon  to  give  values  with  an  absolute  accuracy  approaching  10%.  Beam  foil  spectroscopy 
can  be  included  with  these  more  accurate  methods  if  one  can  be  certain  that  there  is  no 
unrecognized  cascading  from  higher  levels. 

Although  recent  work  has  produced  more  accurate  values  of  some  spectral  lines 
in  selected  elements,  the  work  of  Corliss  and  Bozman  remains  the  broadest  treatment 
of  transition  probabilities  presently  available.  Corliss  and  Bozman  used  the  spectral 
line  intensities  measured  by  Meggers,  Corliss  and  Scribner  to  calculate  transition 
probabilities  of  some  25,000  classified  lines  in  70  elements.  Meggers,  Corliss  and 
Scribner  had  measured  the  peak  intensities,  I,  of  39,000  spectral  lines  radiated  from 
free-burning  10A  dc  arcs  in  air  between  electrodes  of  copper  to  which  1 atom  in  1000 
of  the  clement  of  interest  had  been  added.  Unfortunately  when  the  Corliss  and  Bozman 
values  arc  compared  with  the  more  accurate  values  determined  by  density-independent 
methods,  errors  as  large  as  a factor  of  20  have  been  found.  There  is  therefore  a 
widespread  reluctance  to  use  the  Corliss  and  Bozman  values,  although  in  many  cases 
they  are  the  only  ones  available. 

Because  the  Corliss  and  Bozman  values  remain  the  only  consistent  set  of  transi- 
tion probabilities  spanning  the  stronger  spectral  lines  of  most  of  the  solid  elements, 
there  has  been  considerable  interest  in  attempting  to  improve  these  values.  Three 
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types  of  corrections  have  been  suggested:  an  energy  dependent  correction,  an  inten- 
sity dependent  correction  and  a correction  in  the  copper  arc  temperature.  Huber  and 
Tobey3  and  subsequently  others  have  shown  that  the  energy  dependent  normalization 


function  which  Corliss  and  Bozman  applied  should  in  fact  be  a constant.  Although  orig- 

4 

inally  Corliss  and  Bozman  took  the  copper  arc  temperature  to  be  5100°K,  Corliss 


subsequently  suggested  that  the  temperature  may  actually  be  as  high  as  7000°K.  Others 

5 

have  suggested  temperature  corrections.  Most  recently  Degenkolb  and  Griffiths  sug- 


gested 6000°I<  on  the  basis  of  comparison  with  more  accurate  transition  probabilities 
for  the  three  elements  Mg,  Ca  and  Al. 

In  the  Meggers-Corliss -Scribner  arc  only  one  atom  in  1000  belongs  to  the  ele- 


ment of  interest;  the  rest  are  copper.  When  0.  1 atomic  percent  of  the  copper  atoms 
are  replaced  by  an  element  whose  ionization  potential  (IP)  is  close  to  or  above  that  of 


copper  (IP  = 7.724  eV),  the  physical  properties  of  the  discharge,  such  as  the  electron 
distribution  aud  temperature,  should  not  be  appreciably  changed.  On  the  other  hand  if 
the  ionization  potential  of  the  added  element  is  substantially  lower  than  that  of  copper, 
then  a noticeably  lowering  of  the  electron  and  arc  temperatures  could  be  observed. 

For  an  optically  thin  gas,  one  in  which  effects  due  to  self-absorption  can  be  neg- 
lected, the  intensity  of  the  light  spontaneously  emitted  in  a particular  spectral  line  X 
due  to  N atoms  in  the  upper  energy  level  is 


I = N A 


ui  hc/X 


(0 


If  the  gas  is  in  local  thermodynamic  equilibrium  so  that  the  population  of  the  emitting 
atomic  species  can  be  characterized  by  a Boltzmann  distribution,  then 


N 


8. 
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where  the  g's  represent  the  statistical  weights  of  the  ground  (o)  and  upper  (u)  energy 
levels . 


Now  in  terms  of  the  partition  function 
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and  the  total  number  of  the  emitting  atomic  species 
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we  can  express 


In 


where 


E 

u 


(5) 


n 8tr^e^h  N (6) 

C ' m U 

is  a constant  for  each  light-emitting  species.  Here  the  Ladenburg  oscillator  strength, 
{,  has  been  used  instead  of  the  transition  probability,  Auf.  The  relationship  is 


gf  = g/fu 


8uful 


mcX 

o 2 2 

8n  e 
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If  transition  probabilities  or  f values  are  available  for  an  element  from  one  of 
the  density -independent  methods,  the  presence  of  a straight  line  plot  of  log  (IX  / gf) 
versus  E indicates  that  the  arc  can  be  characterized  by  a Boltzmann  distribution. 
The  arc  temperature  can  then  be  determined  from  the  slope  of  the  line.  When  we 
carried  out  a linear  least  squares  analysis  on  such  plots  for  elements  with  ionization 
potentials  within  1 cV  of  that  of  copper,  we  found  Boltzmann  temperatures  of  approxi 
mately  6200°K  as  indicated  in  Table  I. 

In  fact  we  find  that  if  we  take  account  of  the  different  values  10  for  the  atomic 
partition  functions  by  plotting 


In 


IX_ 

gf 


U = C 


— — E 
kT  u 


(7) 


2 2 

where  C = CU  = - e— - N,  then  the  points  for  all  129  spectral  lines  listed  in  Table  I, 
1 m i ] . 12 

together  with  those  for  2 Ag  spectral  lines  (IP  = 7.  574eV)  and  3 Pb  spectral  lines 

(IP  = 7.  415cV)  can  be  very  satisfactorily  fitted  by  a single  Boltzmann  temperature  of 
6230°K  as  shown  in  Figure  1.  The  fact  that  the  intercept  Cj  is  the  same  for  different 
elements  indicates  that  the  numbers  of  emitting  atoms  are  the  same.  This  corrobo- 
rates Corliss  and  Bozman's  claim  that  the  copper  electrodes  diluted  1 atom  in  1000 
with  the  clement  of  interest  dispense  atoms  of  each  element  into  the  arc  at  the  same 
rate.  Furthermore  it  indicates  that  demixing  effects  which  would  tend  to  reduce  the 
concentration  of  lighter  elements  and  those  with  ionization  potential  substantially  be- 
low that  of  copper  are  not  important  for  the  elements  we  have  examined. 


Contrary  to  the  claim  of  Degcnkolb  and  Griffiths,  Ca  and  A1  do  not  yield  simi- 
lar temperatures.  The  arc  temperature  does  appear  to  decrease  as  the  ionization 
potential  of  the  added  element  in  the  Meggers -Corliss -Scribner  copper  arc  is  lowered 


• Fe  (IP  = 7.87eV) 
XCu  (IP  = 7 724eV) 
-'•Mg  (IP  = 7.644eV) 
x Ag  (IP  = 7.574eV) 
▲ Pb  (IP  = 7.4 1 5 eV) 
■ Ti  (IP  = 6.82  eV) 


ENERGY  OF  UPPER  LEVEL,  Eu  (kK  = IOJcm  ’) 


Fig.  1.  Semilog  plot  of  intensities  of  134  spectral  lines  from  6 elements 
with  ionization  potentials  within  1 eV  of  that  of  copper.  The 
straight  line  least  squares  fit  indicates  a Boltzmann  distribution 
at  6230°K  for  the  Megger s -Cor lis s -Scribner  arc. 
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more  than  approximately  1 eV  below  that  of  copper.  For  example  the  temperature  we 
determined  for  65  spectral  lines  of  Ca  (IP  = 6.  1 1 leV)  was  5037°K  and  for  12  spectral 
lines  of  A1  (IP  = 5.  984eV)  was  544  1°K. 

Therefore  we  suggest  that  for  the  52  elements  in  the  Corliss  and  Bozman  tables 
with  ionization  potentials  > 6.  8eV,  the  temperature  of  the  Meggers -Corliss -Scribner 
arc  should  be  taken  as  6230°K.  In  addition  use  of  the  partition  function  for  each  element 
will  superimpose  the  spectral  lines  of  all  of  these  elements  on  a single  Boltzmann- 
distribution  linear  semilog  plot.  Improved  transition  probabilities  can  then  be  obtained 
without  having  a single  absolute  transition  probability  measured  within  that  element. 
Work  is  continuing  at  the  Polytechnic  to  determine  the  range  of  validity  of  this  proce- 
dure as  well  as  to  determine  the  increase  in  accuracy. 
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COMPUTER  MODELING  OF  THE  DYNAMICS  OF  METAL  VAPOR  LASERS 
N.  Solimene  and  W.  T.  Walter 

The  experimental  investigation  of  metal  vapor  lasers,  in  particular  the  copper 

1 2 
vapor  laser,  and  the  search  for  new  metal  vapor  lasers  are  significant  parts  of  the 

laser  work  at  the  institute.  For  this  reason  it  was  felt  that  numerical  modeling  of  the 
dynamics  of  such  lasers  would  be  a useful  adjunct  to  the  experimental  effort.  Conse- 
quently a start  in  the  development  of  a computer  model  has  been  made. 

A.  The  Model 

The  model  includes  rate  equations  for  the  populations  of  various  energy  levels 
which  are  directly  or  indirectly  involved  in  the  development  of  actual  or  potential  laser 
action.  Also  included  are  rate  equations  for  the  populations  of  additional  atomic  species 
which  may  be  present  to  act  as  buffer  gases  or  to  facilitate  electrical  breakdown  and 
discharge.  The  primary  excitation  is  provided  by  electron  collisions.  De-excitation 
may  bo  by  electron  collisions,  as  well  as  by  other  mechanisms  such  as  radiative  re- 
laxation and  diffusion  to  walls  and  subsequent  de -excitation.  The  rate  constants  for 
transitions  between  the  various  energy  levels  and  for  ionization  depend  on  the  crocs 
sections  for  these  processes  and  on  the  electron  number  density  and  energy  distribu- 
tion. To  simplify  the  model,  the  electron  energy  distribution  is  assumed  to  have  a 
known  functional  form  throughout  the  discharge  pulse,  but  the  average  electron  energy 
is  permitted  to  change  in  a manner  determined  by  an  energy  balance  equation.  Con- 
venient distributions  are  the  Maxwell- Boltzmann  and  the  Druyvesteyn  distribution. 
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This  type  of  model  has  been  used  by  Gerry  to  investigate  the  pulsed  molecular 

4 

nitrogen  laser,  and  by  Leonard  to  investigate  the  copper  vapor  laser. 

Since  in  most  cases  the  cross  sections  for  ionization  or  for  excitation  are  not 

5 

known,  they  have  been  estimated  using  the  semiclassical  formulas  of  Gryzinski. 

The  specific  rate  constants  are  then  calculated  as  a function  of  electron  temperature 
in  the  Boltzmann  case,  or  as  a function  of  the  average  electron  energy  in  the 
Druyvesteyn  case. 

The  circuit  equations  for  capacitor  voltage  and  discharge  current  density  are 
those  used  by  Gerry  and  Leonard.  The  ohmic  resistance  term  depends  on  the  momen- 
tum transfer  cross  sections  for  the  various  species  and  the  electron  energy  distribution. 
It  is  believed  that  the  correct  treatment  of  this  term  will  be  of  some  importance  in 
correlating  experimental  results  involving  different  buffer  species.  This  will  require, 
however,  good  estimates  of  the  momentum  transfer  cross  sections.  For  the  present, 
ad  hoc  estimates  are  being  used.  It  is  hoped  that  classical  phase  shift  calculations 
will  provide  better  estimates  which  will  be  incorporated  later  on. 

The  ohmic  resistance  also  enters  in  an  ohmi  c heating  term  in  the  energy  bal- 
ance equation  for  the  average  electron  energy.  Other  terms  are  energy  transfer  terms 
which  describe  excitation,  de -excitation,  and  ionization.  The  electron  number  density 
is  determined  by  the  ionization  rates  for  the  atomic  species  present  and  by  electron 
loss  terms  such  as  recombination.  The  latter  are  being  ignored  for  the  present  since 
they  are  expected  to  be  important  only  in  the  after-glow  when  the  laser  action  is  com- 
pleted. 

Phenomenological  equations  for  the  laser  photon  densities  are  included  for  each 
of  the  potential  laser  transitions.  Output  power  densities  may  be  computed  using 
effective  optical  cavity  time  constants. 

J3.  The  Equations 

1.  Circuit  Equations 

dV  A . 
dF  = ' C J 

at  = L^(v-ptj )» 

where 


V = voltage  across  capacitor 

j = current  density 

C = capacitance 

L = circuit  inductance 

L - length  of  discharge 

A = cross-sectional  area  of  discharge 

p = discharge  resistivity. 

Note  that  p depends  on  electron  energy  distribution. 
2.  Electron  Equations 
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(k) 


= electron  number  density 
= average  electron  energy 

= atom  number  density  of  component  in  i^1  energy  level 

= ionization  energy  for  k^  species 
= energy  of  i^  level  for  k^  component 

= specific  electron  ionization  rate  for  k^1  component  in  i^1  energy  level 


(k) 


R.V'7  = specific  electron  excitation  (de-excitation)  rate  from  to  level  of 

J . th 

k species, 

(k)  (k) 

Note  that  R.  and  R.^  depend  on  electron  energy  distribution  in  a manner  which  is 
detailed  below. 


3.  Population  Equations 
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k = 1,  2,  c , i=l.  2,  ...,s(k)  ; 


where 


c 

s(k) 

Nik) 

i,  eq 

T(k) 

ij 


number  of  components  in  gas  mixture 

tli 

number  of  energy  levels  for  k component 


= equilibrium  number  density  at  gas  temperature 


relaxation  time  for  i*k  to  level  of  k component  due  to  gas  or 
v/all  collisions 


A?.  = spontaneous  emission  coefficient 

ij 

n(2)  - index  for  upper  level  of  n laser  transition 

^n(2)  = nam^er  density  in  upper  level  of  ntk  laser  transition 

til 

n(l)  = index  for  lower  level  of  n laser  transition 

til  $ ^ 

g.  = degeneracy  of  j energy  level  in  1 component 


6(n,  i) 


+ 1 if  i = n(l ) 

-1  if  i = n{2) 

0 otherwise 


th 


B = stimulated  emission  coefficient  for  n"  laser  transition 
n 

til 

P = photon  density  for  n laser  transition. 


n 

4.  Laser  Photon  Density  Equations 

dP  , / g 

-rP  = - -1  P + B (m(1) 
dt  T n n \ 
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n=  1 1 2,  . . . i n 
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where 


= number  of  laser  transitions 

T = relaxation  time  for  n^1  laser  transition 
n 

V = effective  mode  volume  . 
n 

5.  Specific  Rate  Constants 

Rik>  = f ai(k)<*>  Eo>  dc 


"If  - I <'»  M t(e.  co)  dc  . 


where 


j!k^(e)  = cross  section  for  electron  ionization  of  i1"  energy  level  in  k 
1 component  as  function  of  electron  energy 

o!^(e)  = cross  section  for  electron  excitation  (de-excitation)  from  i to  j 
ij  a.i_ 

energy  level  in  k component  es  function  of  electron  energy 

m = electron  mass 

'(c  c ) ~ electron  energy  distribution  function  with  average  energy  equal  to  Cq. 
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a.  Maxwell- Boltzmann  Distribution 
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= electron  temperatures 


b,  Druyvesteyn  Distribution 


f(c,  r.  = 
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6.  Discharge  Resistivity 
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whe  re : 


e = electron  charge 

(k)  th 

N'  = total  number  density  of  k component 

0<k)(e)  = momentum  transfer  cross  section  for  electron  collisions  with  ktk 
spee’es  as  function  of  electron  energy. 

C.  Numerical  Procedure  and  Results 

The  above  system  of  nonlinear  differential  equations  which  results  from  the 
model  is  integrated  using  the  HPCG  routine  in  the  Fortran  SSP  library. 

The  routine  was,  however,  modified  so  that  the  step  size  is  controlled  by  the  relative 
error  rather  than  absolute  error.  The  reason  for  this  is  the  large  dynamic  range  of 
the  variables  involved  , as  well  as  their  disparate  magnitudes.  Also,  the  allowed 
number  of  step-size  halvings  was  increased  in  order  to  get  past  the  gas-breakdown 
phase  of  the  time  development. 

Some  preliminary  results  have  been  obtained  using  a protolypic  three-level  cop- 
per energy  level  structure  with  no  buffer  species.  These  results  are  reasonable  in 
terms  of  known  experimental  behavior,  but  will  not  be  described,  since  their  only 
purpose  was  to  aid  the  program  development.  It  is  hoped  that  this  model  can  be  suf- 
ficiently well  developed  so  that  it  will  be  capable  of  not  only  correlating  experimental 
data  but  also  of  guiding  future  experiments  in  copper,  bismuth  and  other  metals. 
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APPENDIX  III 
STRATIFIED  WAKES 

The  work  on  the  subject  of  wake  behavior  in  stratified  flows  was 
carried  out  in  prior  contract  periods  (DARPA-STO  program  administered 
by  the  U.S.  Army  Research  Office  under  contract  DAHC04-69-C-007)  . 
However,  the  reports  were  completed  under  the  current  contract  and  will 
be  issued  under  its  auspices. 

One  report  by  G.  H.  Strom^  presents  an  account  of  the  final  develop- 
ment of  the  Polytechnic's  the  rmally-  stratified  air  wind-tunnel,  its  cali- 
bration and  operation.  It  also  presents  experimental  results  of  a series 
of  experiments  on  wakes  in  this  stratified  flow  under  conditions  intended 
to  simulate  wakes  in  ocean  ther moclines , and  related  modeling  criteria. 
This  work  is  believed  to  be  a milescone  achievement  in  stratified  flow 
simulation. 
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A second  report  by  P.  Khosla  , presents  a preliminary  study  of  the 
theoretical  treatment  of  a shear  flow  in  the  form  of  a jet  within  a thermally 
stratified  flow.  It  depicts  some  of  the  difficulties  involved  in  formulating 
a formally  correct  theory,  and  points  the  way  to  further  work  which  may 
lead  to  satisfactory  solutions. 
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